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Abstract 

We propose an exact nonparametric inference scheme for the detection of nonlinear- 
ity. The essential fact utilized in our scheme is that, for a linear stochastic process 
with jointly symmetric innovations, its ordinary least square (OLS) linear predic- 
tion error is symmetric about zero. Based on this viewpoint, a class of linear signed 
rank statistics, e.g. the Wilcoxon signed rank statistic, can be derived with the 
known null distributions from the prediction error. Thus one of the advantages of 
our scheme is that, it can provide exact confidence levels for our null hypothesis 
tests. Furthermore, the exactness is applicable for finite samples. We demonstrate 
the test power of this statistic through several examples. 
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1 Introduction 



Nonlinear statistics such as correlation dimension and Lyapunov exponent 
PQ have been widely adopted in many fields to identify the underlying dy- 
namical systems. However, even for a linear stochastic process with simple 
autocorrelation, these statistics could have finite and predictable values [2J, 
thus without careful treatment, one may mistake a linear stochastic process 
for a nonlinear deterministic system when simply examining whether these 
statistics are convergent. This observation requires us to examine the basic 
properties of the underlying system in order to apply nonlinear analysis meth- 
ods with greater confidence. Based on this viewpoint, to discriminate (station- 
ary) linear stochastic processes from nonlinear deterministic systems, various 
methods have been developed, for example, the idea to investigate the orien- 
tations of the tangents to the system trajectory within given regions [6], the 
proposal to test the continuity of the underlying systems [7] , and the sugges- 
tion to measure the mutual information and redundancy within the framework 
of information theory [8], to name but a few (see [3] for an extensive study). 

In general, these methods focus on exploring the difference of some charac- 
teristic behaviors or properties between linear stochastic and nonlinear (de- 
terministic) systems. In addition, a confidence level is usually preferred in 
order to indicate the reliability of the results. In practical situations, if only 
a scalar time series from an unknown source is available, the conventional 
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approach for inference of the underlying system with the confidence level, as 
proposed by Theiler et al. [4], is to first assign a null hypothesis to the un- 
derlying system, and then apply the bootstrap method to produce a set of 
constrained-realization surrogate data, which should have the same statistic 
distribution as the original time series under the null hypothesis. Hence based 
on the empirical distribution of the adopted statistic of the surrogates and 
the original time series, one then determines whether to reject the null or not. 
However, since the exact knowledge of the statistic distribution is often not 
available, one will resort to certain discriminating criterion to help make the 
decision more objectively and determine the corresponding confidence level (if 
to reject). The popular discriminating criteria appearing in the literature of 
nonlinear (chaotic) time series analysis usually include two classes: parametric 
and nonparametric. 

The parametric criterion assumes that the statistic follows a Gaussian dis- 
tribution, and the distribution parameters, i.e. the mean and the variance, 
are estimated from the finite samples. One can determine whether to reject 
the null by examining whether the statistic of the original time series follows 
the statistic distribution of the surrogates. The corresponding confidence level 
of inference can be calculated from the estimated statistic distribution (see 
our discussions later); The nonparametric criterion [5] examines the ranks of 
the statistic values of the original time series and its surrogates. Supposes 
that the statistic of the original time series is Sq and the surrogate values 
are {sj}^ 1 given N surrogate realizations. Then if the statistic of both the 
original time series and the surrogates follows the same distribution, the prob- 
ability is 1/ (N + 1) for So to be the smallest or largest among all of the values 
{sq, Si...sjv}- Thus if N is large, when one finds that sq is smaller or larger 
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than all of the values in {sj}^, it is likely that so instead follows a dif- 
ferent distribution from that of {sj}^. Hence the criterion rejects the null 
hypothesis whenever the original statistic Sq is the smallest or largest among 
{so, s\...sn}, the false rejection rate is considered as 1/(N + 1) for one-sided 
tests and 2/(N + 1) for two-sided ones. 

Note that the above two simple criteria are often adopted thanks to the dif- 
ficulty in calculating the exact distribution of the nonlinear statistic under 
test. Although these two criteria are heuristic, they are often questionable in 
practice. For example, for the first criterion, the normality assumption may 
not approximate the actual distribution well. For the second, suppose that so 
and {si}^ follow the same distribution and let the range of the surrogate 
statistic be A, while the support of the null distribution be B, and C = B\A 
be the complement set of A given B. Then according to the rule, one rejects 
the null hypothesis whenever the original statistic Sq G C, and the actual false 
rejection rate is the probability Pr(s G Cljsj}^), which usually will not 
simply depend on the number of surrogate realizations. 

In this communication we will propose a new statistic, namely the Wilcoxon 
signed rank statistic, for detection of the potential nonlinearity of a scalar 
time series in the framework of Theiler et al. This statistic could be de- 
rived from the linear prediction error of the time series and proves to be a 
Wilcoxon variate under weak conditions. As an advantage, it has a known 
null distribution, thus inference with exact confidence level becomes possible 
based on the knowledge of the statistic distribution. Furthermore, this statis- 
tic could be applied to a wider range of linear stochastic process, including 
the stationary Gaussian colored noise discussed in [3]. We will introduce the 
detail of the derivation of the Wilcoxon signed rank statistic in Section 2. For 
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demonstration, we will apply this statistic to several examples in Section 3. 
Finally we will summarize the whole work in Section 4. 



2 Methodology 



Now let us begin the introduction of the new statistic. As the first step, we 
would like to specify the null hypothesis Hnto be tested. Given a stationary 
time series {xi} with finite coherence time lj j5j, our objective is to detect if 
there exists nonlinearity of the underlying system. To perform a hypothesis 
test, we instead assume that the time series {x{\ is from a linear stochastic 
process with independent jointly symmetric innovations, which is the null hy- 
pothesis Ho to be tested in our work. For the processes that are not consistent 
with H , they will be attributed to the alternative null hypothesis Hi. For 
clarity, let us explain with more detail (readers are also referred to the origi- 
nal works [lj). In general, the data generation processes can be classified into 
linear or nonlinear stochastic processes and linear or nonlinear deterministic 
systems. However, for linear deterministic systems such as a sine wave, since 
they usually appear very regular (their portraits in phase spaces are fixed 
points or limit cycles), it is not difficult to detect the underlying mechanism. 
Thus in this work we will exclude them from our discussions [f]. 



1 That is, its linear correlation will eventually tend to zero. 

2 A fixed point with disturbance can be considered as a random walk process, which 
is within the scope of our consideration with the constraints imposed in this work. 
For the cases of limit cycles contaminated with noise, their coherence time will 
usually not be finite. For detection of such trajectories, readers are referred to, e.g., 
[20I2T] 
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Also note that, theoretically it is possible that a linear stochastic process has 
asymmetric innovation terms. However, linear stochastic processes with sym- 
metric innovations are often good approximations in many practical situations, 
e.g., the Gaussian distributions of the driving force of a Browning motion, the 
thermal noise in an electronic circuit, the channel noise in telecommunication 
etc. Thus within the scope of our discussion, when we reject the null hypoth- 
esis Hq, we will say that there is a high possibility that the time series is 
generated from a nonlinear system (either deterministic or stochastic) Q. And 
this result can be taken as the primitive step for further investigation of the 
underlying system, for example, model building. 

For any stationary linear stochastic process, it can be presented by an autore- 
gressive moving average (ARMA) process thanks to the Wold's decomposition 
theorem. For convenience in the later discussion, whenever feasible, we will 
use a p-th order autoregressive (AR(p)) processes instead to describe a linear 
stochastic process with the concrete forms of 



where denotes the innovation terms, which are assumed in our null hypoth- 
esis to be independent of Xi, mutually independent of each other and have a 
distribution with joint symmetry. By "joint symmetry" of a stochastic process 
{e^} we mean that there exists some constant \x so that {ej — /i} and {/x — e{\ 
have the same joint distributions, i.e. the probability density function (PDF) 
/(ex — ji, e 2 — n, ...) = /(/i — ei,/i — e 2 , ...) [IT] . Clearly, linear autocorrelated 
Gaussian processes examined in pE] are consistent with our null hypothesis. 

3 Readers are referred to, e.g. \18\ p. 244-245], for the discussion on the formal 
interpretation of the test result. 



Xi = a + V . 

L ✓ 7 
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Yet the coverage could be extended to a wider range, e.g. the linear stationary 
processes with independent (not necessarily identical) and jointly symmetric 
innovations. 

With the above null hypothesis, we then need to choose a discriminating statis- 
tic to determine whether to reject the null or not. To derive our statistic, 
let us first consider the problem of predicting fc-step ahead value of a linear 
stochastic process {x{\ with jointly symmetric innovations {ej}. Let x\ be 
the prediction at time % [x\ = Xi + k if A: ^ 0), and e\ = Xi + k — x\ denote 
the corresponding prediction error. For the purpose of prediction, we choose 
the ordinary least square (OLS) linear predictor x\ = a^o + Y%= 
with a it j being the the (j + l)th coefficient of the OLS predictor, estimated 
based on the history {xj, Xj_2, ...}. In general situations, the OLS predic- 
tor may not perform as well as other procedures like the forward/backward 
least-square algorithm, however, it does possess an interesting property that 
might not be shared by other algorithms. The fact is that, as proved in 
when an OLS linear predictor is used to predict the fc-step ahead value of a 
linear stochastic processes {xi} with jointly symmetric innovations, even if the 
fitting order of the predictor is misspecified (either lower or higher), and thus 
inaccurate estimated parameters are adopted for prediction, the distributions 
of the prediction error e\ will still be symmetric about zero, i.e. e\ and — e\ 
share the same distribution. This fact immediately implies that the probabil- 
ity Pr(ef > 0) = Pr(— e\ > 0) = Pr(ef < 0) = 1/2 when the distribution of e\ 
is continuous so that the probability Pr(ef = 0) = in the sense of Lebesgue 

4 Here p' is the specified fitting order, and {dij : j = 0, 1, are the estimated 

parameters at time i via the criterion of ordinary least squares, which aims to 
minimize the forward squared error of the prediction. 
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measure (see [12] and [T3], Lemma 10.1.24] for more details). 

Let be an indicator series with m data points so that h{e^) = 1 if 

e\ > and ii(e^) = —1 if < 0. Clearly ij is a Bernoulli variate uniformly 
distributed on {-1,1}, i.e. Pr(ii = -1) = Pr(/ i = 1) = 1/2. With this 
knowledge, we can derive a class of linear signed rank statistics 

to test our null hypothesis, where {-S'i(-)} is the set of scores of the series 
{lefl}^ with rank(\ef\) denoting the rankfin the ascending order) of the ab- 
solute value |ef| among {|ef|}™ 1 [I3J p. 252]jj. Here we choose Si(rank(\e^\)) = 
rank(\e!?\) so that 

SRm^YZsX^i) (3) 

is the widely used Wilcoxon signed rank statistic, which is discretely dis- 
tributed. One could obtain the full knowledge of the distribution by enumer- 
ating all of its possible values, however, this is an ineffective way especially 
when m is large. A remedy to this problem is to use the Gaussian distribution 
N(0, m(m + 1) (2m +l)/6) for approximation based on the limit central theo- 
rem. In [Tl"t chaper 2] it is shown that the approximation works well even for 
small numbers, say, m — 6. Therefore in this work we will adopt the strategy 
of distribution approximation. 

Since we know the distribution of the test statistic, based on the realization 
value of SR m , we can determine whether to reject the null hypothesis with 
an exact confidence level. Take two-sided test as an example, if we want the 
type-I error (the false rejection rate of a correct null) to be less than a, then 



5 Since are symmetric about zero, the absolute value of e\ is adopted to remove 
the possible dependence between S{ and I{. 
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we first find two critical values n u and ni such that n u is the largest integer 
satisfying Pi(SR m > n u ) < a/2 and ni is the smallest integer satisfying 
Pr(SR m < rii) < a/2. If for a time series in test, its statistic SR m > n u or 
SR m < rii, then we reject the null hypothesis. The false rejection rate is a, or 
in other words, the confidence level to reject the null hypothesis is 1 — a. The 
procedures to perform one-sided tests are similar, except that we need to locate 
only one critical value, n u or ni, which instead satisfies Pr(S'-R m > n u ) < a or 
Pr (SR m < rii) < a separately for right or left side test. 

For a linear stochastic process consistent with our null hypothesis, in principle 
one would expect that the actual rejection rate will be the same as the nominal 
(pre-specified) one. However, for a nonlinear system violating the hypothesis, 
it would generally engender higher rejection rates because of the asymmetry 
of the prediction errors (especially for the surrogate data, see the discussion 
in the next section), which is the essential idea inside our method. 



3 Numerical results 



We will apply the above idea to test our null hypothesis for nonlinearity de- 
tection. The whole procedures go as follows: For each time series {xi} in test, 
we first predict its one-step ahead values {x}} via the OLS linear predictor 
and calculate the prediction error ^]. Suppose the error series {e t 1 }^ 1 has m 

6 Theoretically any k shall be okay since one always has that Pr(e^ > 0) = 
Pr(— > 0) = Pr(e^ < 0) = 1/2. However, we do recommend the choice of k = 1 
because in practice there might be additional noise from computer (e.g round-off 
effect) based on the OLS predictor if k > 2, which may affect the symmetry of the 
prediction error. 
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data points, then we use the normal distribution N(0, m(m+ l)(2m + 1)/6) to 
find the critical values (n u and ni) for two-sided test at the nominal confidence 
level of 95%. If the calculated Wilcoxon statistic SR m [n h n u ], then we can 
reject the null hypothesis with the false rejection rate of 5%. In the following 
we will demonstrate through several examples the power of our scheme for the 
null hypothesis test. 

The first example is an AR(6) process Xi = S|=i a j x i-j + e i with coeffi- 
cients (ai,...,a 6 ) = (0.6, 0, 0.5, 0, —0.6, 0.3), where innovations {e«} are 
uniformly distributed on interval [0, 0.1] (symmetric about 0.05). The sec- 
ond is an ARMA(1, 1) process, i.e. Xi = a\Xi-\ + e, — &iej_i with parameters 
a\ = b\ — 0.5, where innovation terms {tj} follow the normal distribution 
N(0, 1). The third data generation process (DGP) is the Henon map [TH] 
H(x,y) = (y + 1 — ax 2 , 0.3x), where parameter a is uniformly drawn from 
the interval [1.35, 1.4]. We will take out the first coordinate x for test. The 
final case is the Rossler system [TB] with continuous description equations of 
(x, y, z) = (—y — z, x + 0.15y, 0.2 + xz — cz), where parameter c is uniformly 
drawn from the interval [9.5, 10]. The sampling time is 0.1 time units, and 
the observations for calculation are taken from the second coordinate y. The 
waveforms of the realizations of each DGP are plotted in Fig. [lj 

Since usually one does not know the true order of an underlying process, a 
fitting order has to be specified for prediction. For this purpose, one may adopt 
the Akaike or Schwarz information criterion [22]. However, as aforementioned, 
if the underlying process of the test time series is linear stochastic with jointly 
symmetric innovations, misidentification of the fitting order would also lead 
to the symmetric prediction error. Thus we need not seek the optimal fitting 
order for our prediction. Instead, we could simply choose several fitting orders 
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for all of the DGPs, say, those starting from 6 to 10. To indicate the power 
of our test scheme, for all of the DGPs in examination, we produce 1000 
realizations with 2000 data points for each, and predict the one-step ahead 
values for the last 500 data points. For demonstration, the prediction errors 
with p' = Q are illustrated in Fig. [2] for each DGP (The results with other 
fitting orders are similar and thus not reported here). With the prediction 
errors, we calculate the Wilcoxon signed rank statistic to determine whether 
to reject the null or not. We record the rejection numbers of our null hypothesis 
and indicate them in Table [TJ For the AR(6) and ARMA(1, 1) processes, the 
rejection rates are nearly 5%, as we expect. For the Rossler system, although 
the rejection rates are various for different fitting orders, in general they are 
much higher than 5%, thus we can reject the null hypothesis. While for the 
Henon map, we see that the rejection rates are only slightly higher than the 
nominal one. The explanation may be that, in principle there is no universal 
scores in Eq. (T5]) that bring the most powerful signed rank statistic for all 
systems ([131 Theorem 10.1.19.]), and the sensitivity of this class of statistics 
might dramatically decrease for certain systems. The problem, however, could 
be relieved to some extent when we apply our framework to surrogate tests, 
as to be shown below. 

In practical situations, one often has only a scalar time series on hand. There- 
fore, for the reliability of the test, we suggest that one first uses the bootstrap 
method, such as [I] and [17] , to generate a number of surrogates, and then cal- 
culates the test statistic of the surrogates and determines whether to reject the 
null hypothesis or not. If the actual rejection rate is higher than the nominal 
one (i.e. the rate when the original time series is consistent with our null hy- 
pothesis), then we can safely reject the null hypothesis. For illustration, let us 
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examine the previous examples again. For each example, we generate only one 
sample, and use the bootstrap method to generate 1000 of its surrogates. In 
our test scheme, we adopt the temporal-shift algorithm in [T7f2"T] to generate 
surrogates since this algorithm does not require the Fourier transform and thus 
avoids some of its shortcomings The main idea of the algorithm is that, if a 
time series is linear stochastic with the form of Xt = a + YJj=\ a j x i-j + e %-, 
then for any coefficients (3 and 7, the surrogates {yj = (3xi + 7Xj +r } also 
follow linear stochastic forms with the constants of {(3 + 7)^0 an d the inno- 
vations of + 7ej +r } respectively. A linear stochastic time series {xj} will 
always produce linear stochastic surrogates {yj = (3xi + 7Xj +T }, therefore in 
principle we cannot reject the null hypothesis test via the Wilcoxon statistic. 
Following the suggestion in [2T], in our calculations we let parameter j3 be 
uniformly drawn from the interval [0.6, 0.8] and parameter 7 = (1 — 2 ) 1 / 2 to 
produce surrogates. The results of null hypothesis test are presented in Table 
El which appear different from those in Table [U For the two linear stochas- 
tic processes, the rejection numbers decrease to zero. This is because that 
the original samples of the linear stochastic processes are not rejected in the 
tests (as 95% of the samples are), therefore all its surrogates, which are es- 
sentially the addition of two segments of the original samples, will also not be 
rejected in the tests. However, the situation for nonlinear systems is different. 
Adding together two segments of the original samples from nonlinear systems 
increases the complexity, thus the rejection rates of the surrogates increase 
correspondingly compared to those in Table [TJ 

For comparison, we also generate 1000 surrogates through the constrained- 



7 A good review of the conventional surrogate algorithms can be found in, for 
example, [181 chapter 11] 
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realization method [I], which preserves the linear correlation of the original 
time series but introduce some randomization to the phase of the Fourier 
transform. The computation results are indicated in Table [3J from which we 
could see that, the rejection rates of linear cases drop as well in contrast to 
those in Table [TJ However, the rates of nonlinear cases remain close to those 
in Table [U this is because, with the preservation of the linear correlation, 
the phase perturbation of the surrogates does not significantly affect the data 
structure, thus the rejection rates of the nonlinear systems are close to those 
in Table[TJ 

An further question is that: what will happen if the data in test is linear 
but with asymmetric innovation terms? To answer this question, let us see an 
example. We construct the following process 



with the innovation term e« following the beta(2, 5) distribution. Here 



with B(a,/3) denoting the beta function. We choose the same computation 
settings as previously adopted. The results are summarized as follows, 

(1) For 1000 realizations generated by Eq. (j3J), the rejection numbers corre- 
sponding to the fitting orders (6 to 10) are 98, 99, 103, 105, 109 respec- 
tively. 

(2) For 1000 surrogates produced by the temporal shift algorithm based on 
one realization of Eq. (jllh the rejection numbers corresponding to fitting 
orders are 0, 0, 0, 0, respectively. 

(3) For 1000 surrogates produced by the conventional constrained realization 



Xi = 0.1 + 0.8xj_i — 0.5xj„ 2 + 0.2xj_ 3 + ej 



(4) 
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algorithm, the rejection numbers are 109, 97, 100, 117, 129 respectively. 

Based on the results, one can see that in (3) one can correctly reject the null 
hypothesis but in (2) one fails. Because the temporal shift algorithm is simply 
based on the linear superposition principle, although failing to reject the null 
hypothesis in (2), it does tell us some useful information, i.e., the data in 
test is unlikely to be nonlinear (otherwise the rejection rate might obviously 
exceed the nominal, 5%, because of the superposition of the surrogates). In 
contrast, rejecting the null hypothesis in (3) says that the data in test cannot 
be generated by a linear process with jointly symmetric innovation terms. 
Thus conbining the results in (2) and (3) might shed light on the underlying 
process, which is to be investigated in a future work. 

As an application, we apply our scheme to examine an electrocardiogram 
(ECG) record during ventricular fibrillation (VF), a segment of which is plot- 
ted in Fig. [3] (for the details of data acquisition, see [19]). The segment in 
examination has 50, 000 data points, thus we use the temporal-shift method to 
generate 1000 surrogates since it needs less computations than the constrained- 
realization algorithm. We also adopt the same fitting orders, from 6 to 10, to 
predict 100 one-step ahead values. As the results, the correspondingly rejec- 
tions of our null hypothesis are 566, 541, 530, 550, and 521 (all out of 1000 
tests), a strong hint of nonlinearity. 

4 Conclusion 

To summarize, we have proposed an exact nonparametric inference scheme 
to detect the potential nonlinearity in a scalar time series. The exactness 
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of our inference comes from the knowledge of the exact distribution of the 
adopted discriminating statistic, i.e., the Wilcoxon signed rank statistic, which 
indicates remarkable test power through the several examples examined in 
this communication. The advantages of this statistic include: it possesses the 
known null distribution. Thus it is easy for us to find the exact confidence 
interval for the inference of the underlying process. Furthermore, the exactness 
of the statistic distribution does not rely on the size of the sample time series 
in test. Comparatively, for many nonlinear discriminating statistics adopted 
in the literature, e.g. the correlation dimension, computations with too short 
samples may cause serious distortions. 

XL was supported by a Hong Kong University Grants Council Competitive 
Earmarked Research Grant (CERG) number PolyU 5216/04E. 
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a) One-step ahead prediction errors of the AR(6) realization 




(b) One-step ahead prediction errors of the ARMA(1 ,1 ) realization 



100 200 300 400 500 




100 200 300 400 500 



(c) One-step ahead prediction errors of the Henon map realization 
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(d) One-step ahead prediction errors of the Rossler system realization 
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Fig. 2. One-step ahead prediction errors (p f = 6) of a realization of: (a) the AR(6) 
process, (b) the ARMA(1, 1) process, (c) the Henon map, and (d) the Rossler 
system. 
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Fig. 3. Waveform of a human electrocardiogram (ECG) record during ventricular 
fibrillation (VF) 
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Table 1 

Numbers of rejections of the null hypothesis (out of 1000 
replica) for data generation processes (DGPs) with different 
fitting orders. 

DGP Rejections for the fitting order of 

6 7 8 9 10 
AR(6) 55 54 54 53 50 
ARMA(1,1) 44 41 43 44 48 
Henon 76 55 60 62 67 
RSssler 990 442 752 298 458 
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Table 2 

Rejections of the null hypothesis for 1000 temporal-shift sur- 
rogates of DGPs with different fitting orders. 

DGP Rejections for the fitting order of 

6 7 8 9 10 
AR(6) 

ARMA(1,1) 
Henon 432 307 308 203 207 
RSssler 1000 553 1000 593 626 
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Table 3 

Rejections of the null hypothesis for 1000 constrained- 
realization surrogates of DGPs with different fitting orders. 

DGP Rejections for the fitting order of 

6 7 8 9 10 
AR(6) 14 14 10 11 9 
ARMA(1,1) 112 1 2 
Henon 83 56 63 61 83 
RSssler 992 449 753 301 460 
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